Анализ номинативных данных

Проверка гипотезы о распределении номинативной переменной

1. Хотим проверить гипотезу о том, что распределение некоторой номинативной переменной отличается от заданной.

Предположим, какая партия является фаворитом.

Гипотеза — наше империческое распределение частот, отличается от заданного теоретического.

Нулевая гипотеза — все три партии встречаются равновероятно. Однако очевидно мы видим что в нашей выборке существуют эмпирические отклонения.

Первое что научимся делать — подсчитывать на сколько значимым является это отклонение

Проверка гипотезы о распределении номинативной переменной

Проверка гипотезы о взаимосвязи двух номинативных переменных

Дaнные о том, подействовало ли лекарство на больного

Проверка гипотезы о взаимосвязи двух номинативных переменных

Обе переменные номинативные — лекарство а, лекарство б. Поправился или не поправился

Более сложные модели

Хотим предсказать, как влияют другие переменные на выздаровление.

Зависимая переменная нашего анализа будет являться номинативная переменная с двумя градациями. Такой тип регрессии называется логистическая регрессия.

Подброс монетки

Результат бросания монетки это — номинативная переменная орёл или решка. Подбрасываем 60 раз. Предположим орёл выпал 40 раз, решка 20.

Нулевая гипотеза: распределение орла и решки равновероятна.

Альтернативная гипотеза: распределение орла отлично от распределения решки. Иными словами

\[H_0: p_{орла} = 0.5\] \[H_1: p_{орла} \neq 0.5\]

Первое что мы должны сделать — проверить на сколько отклоняются предсказанные от полученных.

Наблюдаемые обозначим буквой O (англ. Observe). Предсказанные E (Expected).

Итак как работает статистика — она говорит, окей есть нулевая гипотеза, есть ожидаемые резуьтаты, есть эмпирические данные. Смотрим есть ли расхождения. Нужно получить некое численное значение.

Далее мы говорим: а если бы нулевая гипотеза была верна, и мы много раз повторяли эксперимент, то как бы вела некая интересующая нас статистика.

Первое что приходит на ум: смотреть разность ожидаемых и наблюдаемых.

\[d = (O_1-E_1) + (O_2-E_2) = (40-30) + (20-30) = 0\]

Это не самый лучший вариант, посколько в итоге мы получили что отклонений нет

Можно возвести в квадрат разности, чтобы избавиться от минуса.

\[d^2 = (O_1-E_1)^2 + (O_2-E_2)^2 = (40-30)^2 + (20-30)^2 = 200\] Но такой вариант тоже плохой. Потому что в такой формуле получается, что одно и то же значение \(d^2\) может характеризовать как довольно значительные, так и практически несущественные расхождения между наблюдаемыми и ожидаемыми частотами. Если бы у нас была другая выборка:

То мы получили бы точно такое же численное значение

\[d^2 = (O_1-E_1)^2 + (O_2-E_2)^2 = (1040-1030)^2 + (1020-1030)^2 = 200\]

Пирсон тоже понял что такой подход неправильный. И математически доказал, что корректна будет формула

\[X^2_{Pearson} = \sum_{i=1}^{k}\bigg(\frac{O_i-E_i}{\sqrt{E_i}}\bigg)^2\] k — это все наши ячейки. Из этой формулы следует два простых вывода:

• Расстояние хи-квадрат Пирсона может равняться нулю.

• Расстояние хи-квадрат не может быть отрицательным.

Или в нашем случае это будет

\[X^2_{Pearson} = \sum_{i=1}^{k}\bigg(\frac{O_i-E_i}{\sqrt{E_i}}\bigg)^2 = (\frac{(40-30)}{\sqrt{30}})^ 2 - (\frac{(20-30)}{\sqrt{30}})^2 = \frac{100}{30} + \frac{100}{30} \approx 6.7 \]

Теперь второй вопрос — если верна нулевая гипотеза, как бы вёл себя этот показатель \(\chi^2\) при многократном повторении. И для этого нам надо познакомиться с распределением \(\chi^2\)

Пример

В одном из опытов эмпирическое распределение частот некоторого цвета гороха приняло следующий вид: \(O = \{18,55,27\}\)

Чему будут равны ожидаемые значения частот, если предполагаемое теоретическое распределение имеет следующий вид: 1:2:1

Решение: У нас есть суммарное количество появления разных признаков \(18+55+27=100\). Есть суммарное отношение ожидаемых частот \(1+2+1=4\). Получаем что одна часть ожидаемого признака равняется \(25\). Значит ожидаемые значения частот равны \(E = \{25,50,25\}\)

Пример 2

Рассчитаем для этого случая значение \(\chi^2\)

\[X^2_{Pearson} = (\frac{(18-25)}{\sqrt{25}})^2 + (\frac{(55-50)}{\sqrt{50}})^2 + (\frac{(27-25)}{\sqrt{25}})^2 = \frac{196+50+16}{100}=2.62\]

Симуляция распределения

Итак что будем симулировать:

1. 60 раз подбрасываем монетку.

2. Записываем сколько выпало орлов, а сколько решек.

3. Ожидаем что выпадет 30 на 30.

4. Считаем критерий хи-квадрат.

5. Повторяем пункты 1-4 100 раз

100 раз повторяется эксперимент с подбрасыванием симметричной монетки 60 раз. Расстояние хи-квадрат рассчитывается для каждого эксперимента:

chi_square_distance <- function() {
  e <- c(30, 30)
  sum((table(sample(x = 2, size = 60, replace = TRUE)) - e) ^ 2 / e)
}

results <- replicate(100, chi_square_distance())
hist(results)

Так замечательно, а теперь мы хотим знать — какова вероятность получить такие или ещё более выраженные отклонения.

Если верна нулевая гипотеза — то чаще всего отклонения между ожидаемым числом и наблюдаемым не будут большие. Сильные отклонения будут встречаться реже.

Это говорит, о том, что распределение будет нормальным

Распределение Хи-квадрат

Определение распределения Хи-квадрат

Мы вяснили — если мы многократно повторяем эксперимент, и смотрим не сумма, а каждое слагаемое в отдельности. То оказывается, что если верна нулевая гипотеза, то отклонения чаще всего не значительны и будут встречаться равновероятно как в одну так и в другую сторону.

Распределение Хи-квадрат с k-степенями свободы — это распределение суммы квадратов k-независимых стандартных нормальных случайных величин.

Стандартное распределение — нормальное распределение, где среднее равно нулю, а дисперсия единице.

Независимые — это тот факт, что каждое последующее измерение не зависит от измерения предыдущего.

Для того чтобы понять откуда берётся каждое значение в распределении Хи-квадрат с двумя степенями свободы, давайте визуализируем весь процесс.

Предположим мы сделали программу, которая генерирует две стандартные случайные величины.

Вероятность получить значение хи-квадрат равное или большее 5.9 равняется 0.05. Это и есть критическое значение для

Распределение Хи-квадрат

При достаточно большом количестве степеней свободы, распределение \(\chi^2\) стремится к нормальному.

Распределение Хи-квадрат

Задача

Какой процент наблюдений лежит в диапазоне от 2 до 4 у распределения хи-квадрат с двумя степенями свободы?

Калькулятор p-уровня значимости

Ответ:

• Для 4 — P(X > 3.99) =

• Для 2 — P(X > 2) = 0.368

Ответ: \[0.368 - 0.136 = 0.233\]

Задача

Теперь рассчитаем p-уровень значимости для нашего примера с игральной костью. Напомню, что мы получили следующие значения наблюдаемых частот (от единички до шестерки):

10,10,10,5,10,15

Проверьте нулевую гипотезу о том, что эмпирическое распределение частот не отличается от равномерного. В поле для ответа введите получившийся p-уровень значимости.

Калькулятор p-уровня значимости

Решение * df = 5

• O = \(\{10,10,10,5,10,15\}\)

• E = \(\{10,10,10,10,10,10\}\)

\[\chi^2 = 4 \times \bigg(\frac{10-10}{\sqrt{10}}\bigg)^2 + \bigg(\frac{5-10}{\sqrt{10}}\bigg)^2 + \bigg(\frac{15-10}{\sqrt{10}}\bigg)^2 = 0 + \frac{25}{10} +\frac{25}{10} = 5\]

Ответ

Задача

Вернемся к нашему примеру с политическими партиями! Проверьте гипотезу о том, что в ГС нет никаких различий в предпочтениях трех партий. Введите в поле для ответа получившееся значение статистики хи-квадрат с точностью хотя бы до одной цифры после запятой.

• партия А = 10 голосов

• партия Б = 30 г.

• партия В = 50 г.

df = 2

Нулевая гипотеза — распределение для всех партий одинаково. Т.е. если всего 90 голосов, то

• E = \(\{30, 30, 30\}\)

• O = \(\{10, 30, 50\}\)

\[\chi^2 = \bigg(\frac{10-30}{\sqrt{30}}\bigg)^2 + 0 + \bigg(\frac{50-30}{\sqrt{30}}\bigg)^2 =\frac{80}{3} = 26.67\]

Проинтерпритиуем полученный рузельтат.

p < 0.05 — отклоняем нулевую гипотезу.

Принимаем альтернативную гипотезу, что распределение предпочтений избирателей отличается от равномерного. Отвергаем нулевую гипотезу о том, что число сторонников каждой из трех партий в генеральной совокупности одинаково.

(0.53 * 1500) - 1500

## [1] -705

45 * 45 * 2 / 750

## [1] 5.4

Задача

В 2013 году Эдвард Сноуден передал СМИ секретную информацию АНБ, касающуюся слежки американских спецслужб за информационными коммуникациями между гражданами. Однако его поступок вызвал неоднозначную реакцию в обществе. Исследовательский центр USA TODAY провел опрос 1500 граждан США с целью выяснить, воспринимают ли они поступок Сноудена как положительный или отрицательный. 53% опрошенных респондентов оценили разоблачение положительно.

При помощи теста хи-квадрат проверьте нулевую гипотезу о том, что в генеральной совокупности распределение отношения к поступку Сноудена является равномерным, то есть 50 на 50.

Решение

• n = 1500

• df = 1

• E = \(\{750; 750\}\)

• O = $ {0.53 * 1500; 0.47 * 1500} = {795; 705}$

• Нулевая гипотеза — в генеральной совокупности распределение отношения к поступку Сноудена является равномерным, то есть 50 на 50.

• Альтернативная — распределение отношения является не равномерным

\[\chi^2 = \bigg(\frac{795-750}{\sqrt{750}}\bigg)^2 + 0 + \bigg(\frac{705-750}{\sqrt{750}}\bigg)^2 =\frac{4050}{750} = 5.4\]

• \(P(X > 5.4) = 0.0201\) — отвергаем нулевую гипотезу, принимаем альтернативную.

Расчёт ожидаемых значений

Как тогда заполнить таблицу ожидаемых значений? Можно просто поделить общее число человек на 4 и в каждой ячейке записать по примерно 10 человек. Но это плохой подход.

Тот факт что две переменные не зависят между собой, надо интерпритировать следующим образом — распределение соотношения мужчин и женщин как для биологов так и для информатиков.

Например если распределение Ж-М для информатиков 30 на 70. То и для биологов должно быть такое же распределение 30 на 70.

Если у нас дисбаланс классов, мы должны заполнить табличку сопряжённости так, как если бы распределение признаков было в классах одинаково.

Расчет ожидаемых значений

Итак, теперь в общем виде это выглядит так:

Расчет ожидаемых значений

Где:

• \(f_i\) — число наблюдений в \(i\)−ой строке

• \(f_j\) — число наблюдений в \(j\)−ом столбце

• \(N\) — общее число наблюдений в таблице

Переведём это в программные вычисления

O <- matrix(c(15, 11, 9, 6), ncol = 2)
E <- outer(rowSums(O), colSums(O)) / sum(O)
# выводим результаты
print(addmargins(O))

##           Sum
##     15  9  24
##     11  6  17
## Sum 26 15  41

print(addmargins(E))

##                        Sum
##     15.21951  8.780488  24
##     10.78049  6.219512  17
## Sum 26.00000 15.000000  41

# готовое решение
print(chisq.test(O)$expected)

##          [,1]     [,2]
## [1,] 15.21951 8.780488
## [2,] 10.78049 6.219512

# Задача
O1 <- matrix(c(10, 5, 6, 15), ncol = 2)

chisq.test(O1)$expected

##          [,1]      [,2]
## [1,] 6.666667  9.333333
## [2,] 8.333333 11.666667

Расчёт значения Хи-квадрат

Само вычсиление не должно вызывать трудности.

Расчет ожидаемых значений

Основной вопрос в том — сколько степеней свободы у этого распределения.

С одной стороны можно подумать, что здесь 4 степень свободы — по числу ячеек. Но это на самом деле не верно.

Потому что зная только одно значение ячейки и сумму во всех остальных ячейках мы можем вычислить значение в каждой ячейке.

Это означает что корректной t-статистикой будет значение степеней свободы равное 1.

\[df = 1\]

Поправка Йетса

Расчет ожидаемых значений

Расчёт P-уровня значимости

На слайде ошибка, но ход мыслей один и тот же.

Расчет ожидаемых значений

Вот корректное вычисление:

O <- matrix(c(15, 11, 9, 6), ncol = 2)
E <- outer(rowSums(O), colSums(O)) / sum(O)
test <- sum((O - E)^2 / E)
df <- (ncol(O) - 1) * (nrow(O) - 1)
pval <- pchisq(test, df, lower.tail = FALSE)
# сравним результат
print(list(chisq = test, df = df, p.value = pval))

## $chisq
## [1] 0.02087104
## 
## $df
## [1] 1
## 
## $p.value
## [1] 0.8851308

print(chisq.test(O, correct = FALSE))

## 
##  Pearson's Chi-squared test
## 
## data:  O
## X-squared = 0.020871, df = 1, p-value = 0.8851

Итоги

Критерий Хи-квадрат Пирсона это метод который позволяет нам анализировать номинативные данные.

Первый тип гипотезы — распределение частот некоторого признака отличается от заданного теоретического.

Второй — две переменные взаимосвязаны между собой.

Расчет ожидаемых значений

Этот критерий говорит о поведении данных в сумме, но мы не можем сказать в какой именно ячейке было отклонение.

Например в случае голосования за партии — единственное что мы можем сказать, что распределение выбора партий отличается от равномерного.

Задача

Перед знаком стоп некоторые водители останавливаются полностью, другие лишь сбавляют скорость, но некоторые не останавливаются вообще. Важнейший вопрос, есть ли взаимосвязь между полом и стилем вождения автомобиля! Ниже представлена таблица сопряженности данных исследования, посвященного этому вопросу.

O <- matrix(c(20, 11, 7, 15, 12, 9), ncol = 2)
E <- outer(rowSums(O), colSums(O)) / sum(O)
test <- sum((O - E)^2 / E)

print(chisq.test(O, correct = FALSE))

## 
##  Pearson's Chi-squared test
## 
## data:  O
## X-squared = 0.95441, df = 2, p-value = 0.6205

Задача

Обратимся к данным о катастрофе «Титаника». На графике представлена взаимосвязь пола пассажира «Титаника» и того, выжил он или нет в катастрофе. Размеры прямоугольников отвечают за пропорции наблюдений в той или иной ячейке, а цвет прямоугольников — за значение стандартизованного остатка в ячейках. Какие выводы мы можем сделать, проанализировав данный график?

Расчет ожидаемых значений

Ответ:

• Есть все основания отклонить нулевую гипотезу об отсутствии взаимосвязи пола и вероятности выжить в катастрофе.

• Значимые отклонения между наблюдаемыми и ожидаемыми результатами получены во всех ячейках, что позволяет говорить: мужчины вероятнее погибнут, чем выживут, а женщины — наоборот.

• На борту Титаника большинство пассажиров - мужчины

Задача

А теперь изучим данные о взаимосвязи шанса выжить в кораблекрушении «Титаника» и класса билета пассажира. В данном случае график стандартизированных остатков построен по результатам таблицы сопряженности 3 на 2.

Какие выводы мы можем сделать в данном случае?

Расчет ожидаемых значений

Ответ:

• Для пассажира с билетом первого класса вероятность выжить выше вероятности погибнуть.

• Наши данные не позволяют нам сделать вывод о статистически значимом различии в вероятности выжить или погибнуть у пассажиров второго класса.

• Пассажиры из третьего класса чаще погибали в катастрофе, чем пассажиры второго класса.

Задача

Допустим, мы решили провести исследование, целью которого было доказать влияние интересов водителей на безопасность их вождения. На выборке из 100 человек мы зафиксировали следующие показатели: смотрит ли человек онлайн-курсы на Stepic и попадал ли водитель в ДТП за последний год. В результате мы получили следующие данные:

#воссоздадим таблицу
dtp <- rbind(c(10, 40), c(35, 15))
#подпишем строки и столбцы
colnames(dtp) <- c("Был в ДТП", "Не был в ДТП")
rownames(dtp) <- c("Проходит курсы", "Не проходит курсы")
#вот наш график
mosaicplot(dtp, color=T, shade=T, xlab="ДТП", ylab="Group", cex.axis=1, main="")

test <- chisq.test(dtp)
test

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  dtp
## X-squared = 23.273, df = 1, p-value = 1.406e-06

Ответ

Не стоит забывать, что ошибка корреляции может быть совершена и при использовании критерия хи-квадрат. Полученные данные говорят лишь о взаимосвязи двух переменных, но при таком дизайне исследования мы не можем быть уверены в наличии причинно-следственной связи между переменными!

Как работает точный критерий Фишера

Использем запись в общем виде

Расчет ожидаемых значений

Нулевая гипотеза — вероятность поправиться от лекарства 1 = вероятности лекарства 2. Нет никаких различий между использованием двух лекарств. Допустим эта вероятность = 0.5

Тогда можно ввести следующее обоззначение

X — кол-во положительных исходов у пациентов №1

Y — кол-во положительных исходов у пациентов №2

Какова вероятность, что мы веберем 10 пациентов (Л1), мы знаем что вероятность полож исхода = 0.5. Какова вероятность, что мы будем набюлюдать троих человек с положительным исходом.

Это биноминальное распределение. Уравнение Бернулли нам поможет.

Если мы хотим ответить на поставленный вопрос. Вероятность будет равняться по формулле Бернулли, на скрине выше

Итак в общем виде это будет выглядеть вот так

Расчет ожидаемых значений

Теперь если говорить про наш пример, то всё будет выглядеть вот так:

Для вычисления биномиальных коэффициентов (сочетаний) в R используется функция choose (n, k)

choose(4, 3)*choose(4, 1)/choose(8, 4)

## [1] 0.2285714

Итоги

• Критерий хи-квадрат — Анализ таблицы сопряженности произвольного размера, где значения ожидаемых и наблюдаемых частот в каждой ячейке больше 10

• Критерий хи-квадрат с поправкой Йетса — Анализ таблицы сопряженности два на два, где значения ожидаемых и наблюдаемых частот в каждой ячейке больше 5, но меньше 10

• Точный критерий Фишера — Анализ таблицы сопряженности, в которой экстремально мало наблюдений, в некоторых ячейках число ожидаемых и наблюдаемых частот может быть меньше пяти

Задача

Напишите функцию smart_test, которая получает на вход dataframe с двумя номинативными переменными с произвольным числом градаций. Функция должна проверять гипотезу о независимости этих двух переменных при помощи критерия хи - квадрат или точного критерия Фишера.

Если хотя бы в одной ячейке таблицы сопряженности двух переменных меньше 5 наблюдений, функция должна рассчитывать точный критерий Фишера и возвращать вектор из одного элемента: получившегося p - уровня значимости.

Если наблюдений достаточно для расчета хи-квадрат (во всех ячейках больше либо равно 5 наблюдений), тогда функция должна применять критерий хи-квадрат и возвращать вектор из трех элементов: значение хи-квадрат, число степеней свободы, p-уровня значимости.

smart_test <-  function(x){
  y <- table(x)
test <- sum(sapply(y >= 5, sum)) == length(sapply(y >= 5, sum))
ifelse(test, 
       return(c(chisq.test(y)$statistic, chisq.test(y)$parameter, chisq.test(y)$p.value)),
       return(fisher.test(y)$p.value))
}

Итак, проверим функцию

test_data <- as.data.frame(list(am = c(1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), vs = c(0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1)))

smart_test(test_data)

## X-squared        df           
## 0.3475355 1.0000000 0.5555115

Мне понравилось ещё вот такое решение:

smart_test <-  function(x){
  t <- table(x)
  ifelse(any(t < 5),
   return(fisher.test(t)$p.value),
   return(unlist(chisq.test(t)[1:3]))
   )
}

Задача

Почувствуй себя биоинформатиком! Вся наследственная информация в живых организмах хранится внутри молекулы ДНК. Эта молекула состоит из последовательности четырех “букв” — A, T, G и C.

Напишите функцию most_significant, которая получает на вход dataframe с произвольным количеством переменных, где каждая переменная это нуклеотидная последовательность.

В этом примере три последовательности V1 , V2, V3. Для каждой переменной мы можем проверить нулевую гипотезу о том, что все нуклеотиды (A, T, G, C) встречаются равновероятно внутри этой последовательности. Однако, возможно, что в некоторых последовательностях распределение частоты встречаемости каждого нуклеотида отличается от равномерного.

Функция должна возвращать вектор с названием переменной (или переменных), в которой был получен минимальный p - уровень значимости при проверке гипотезы о равномерном распределении нуклеотидов при помощи критерия хи - квадрат.

DNA_generator <- function(n, ncol){
  data <- replicate(ncol,sample(c("A","T","G","C"), n, replace = T))
  rownames(data) <- 1:nrow(data)
  colnames(data) <- paste0("V", 1:ncol)
  return(data)
}

test_data <- as_tibble(DNA_generator(30, 3))

most_significant <-  function(x){
variable_p_values <- sapply(test_data, 
                            function(x) chisq.test(table(x))$p.value)
name_of_minimal <- colnames(test_data[which(variable_p_values == min(variable_p_values))])
return(name_of_minimal)
}

most_significant(test_data)

## [1] "V3"

most_significant <- function(df) {
  result <- sapply(df, function(x) chisq.test(table(x))$p.value)
  names(result)[result == min(result)]
}

Задача

Создайте новую переменную important_cases - фактор с двумя градациями (“No” и “Yes”). Переменная должна принимать значение Yes, если для данного цветка значения хотя бы трех количественных переменных выше среднего. В противном случае переменная important_cases будет принимать значение No.

vector <- c(5,9)
matrix_test <- matrix(c(4,10,4,13,2,8),3,2)

t(t(matrix_test) < vector)

##       [,1]  [,2]
## [1,]  TRUE FALSE
## [2,] FALSE  TRUE
## [3,]  TRUE  TRUE

# Получаем вектор средних количественных значений
mean_iris <- sapply(iris[1:4], mean)

# Сравниваем матрицу построчно с вектором средних
compare_matrix <- t(t(iris[1:4]) > mean_iris)

# Получаем вектор с количеством измерений выше средних
count_high_mean <- rowSums(compare_matrix)

# Задаём условие --- если 3 или больше измерений больше среднего присвой "Yes", иначе "No"

important <- ifelse(count_high_mean >= 3, "Yes", "No")

# Делаем этот вектор фактором

as.factor(important)

##   [1] No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  No 
##  [18] No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  No 
##  [35] No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  No  Yes
##  [52] Yes Yes No  Yes No  Yes No  Yes No  No  Yes No  Yes No  Yes No  No 
##  [69] Yes No  Yes Yes Yes Yes Yes Yes Yes Yes Yes No  No  No  No  Yes No 
##  [86] Yes Yes Yes No  No  No  Yes No  No  No  No  No  Yes No  No  Yes No 
## [103] Yes Yes Yes Yes No  Yes Yes Yes Yes Yes Yes No  No  Yes Yes Yes Yes
## [120] Yes Yes No  Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes
## [137] Yes Yes Yes Yes Yes Yes No  Yes Yes Yes Yes Yes Yes Yes
## Levels: No Yes

# Присваиваем в новую переменную в iris

iris$important_cases <- as.factor(important)

iris$important_cases <- as.factor(ifelse(rowSums(t(t(iris[1:4]) > sapply(iris[1:4], mean))) >= 3, "Yes", "No"))

str(iris$important_cases)

##  Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...

table(iris$important_cases)

## 
##  No Yes 
##  81  69

Задача

Обобщим предыдущую задачу! Напишем функцию get_important_cases, которая принимает на вход dataframe с произвольным числом количественных переменных (гарантируется хотя бы две переменные). Функция должна возвращать dataframe с новой переменной - фактором important_cases.

Переменная important_cases принимает значение Yes, если для данного наблюдения больше половины количественных переменных имеют значения больше среднего. В противном случае переменная important_cases принимает значение No.

Переменная important_cases - фактор с двумя уровнями 0 - “No”, 1 - “Yes”. То есть даже если в каком-то из тестов все наблюдения получили значения “No”, фактор должен иметь две градации.

test_data <- data.frame(V1 = c(16, 21, 18), 
V2 = c(17, 7, 16), 
V3 = c(25, 23, 27), 
V4 = c(20, 22, 18), 
V5 = c(16, 17, 19))



get_important_cases <- function(df){
df$important_cases <- ifelse(rowSums(t(t(df) > sapply(df, mean))) > length(df)/2, 1, 0)
df$important_cases <- factor(df$important_cases, levels = c(1, 0), labels = c("Yes", "No"))
return(df)
}

get_important_cases(test_data)

##   V1 V2 V3 V4 V5 important_cases
## 1 16 17 25 20 16              No
## 2 21  7 23 22 17              No
## 3 18 16 27 18 19             Yes

Задача

Напишите функцию stat_mode, которая получает на вход вектор из чисел произвольной длины и возвращает числовой вектор с наиболее часто встречаемым значением. Если наиболее часто встречаемых значений несколько, функция должна возвращать несколько значений моды в виде числового вектора.

v <- c(1, 2, 3, 3, 3, 4, 5)
v1 <- c(1, 1, 1, 2, 3, 3, 3)
v2 <- c(5,9,20,8,11,7,2,7,7,18,17,15,2,11)


stat_mode <- function(v){
  result <- names(which(table(v) == max(table(v))))
  return(as.integer(result))
  }

stat_mode(v2)

## [1] 7

Задача

Доктор Пилюлькин решил вооружиться статистикой, чтобы сравнить эффективность трех лекарств! Давайте поможем ему и напишем функцию max_resid, которая получает на вход dataframe с двумя переменными: типом лекарства и результатом его применения.

Drugs - фактор с тремя градациями: drug_1, drug_2, drug_3.

Result - фактор с двумя градациями: positive, negative.

Функция должна находить ячейку таблицы сопряженности с максимальным значением стандартизированного остатка и возвращать вектор из двух элементов: название строчки и столбца этой ячейки.

test_data <- read.csv("https://stepic.org/media/attachments/course/524/test_drugs.csv")

max_resid <- function(x){
  table_resid <- chisq.test(table(x))$stdres
  result <- c(rownames(which(table_resid == max(table_resid), arr.ind = TRUE)), 
    rownames(which(t(table_resid) == max(table_resid), arr.ind = TRUE)))
  return(result)
}

max_resid(test_data)

## [1] "drug_1"   "positive"

Задача

library(tidyverse)
str(diamonds)

## Classes 'tbl_df', 'tbl' and 'data.frame':    53940 obs. of  10 variables:
##  $ carat  : num  0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 2 6 7 7 6 5 2 5 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 2 3 5 4 2 6 7 3 4 5 ...
##  $ depth  : num  61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
##  $ table  : num  55 61 65 58 58 57 57 55 61 61 ...
##  $ price  : int  326 326 327 334 335 336 336 337 337 338 ...
##  $ x      : num  3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
##  $ y      : num  3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
##  $ z      : num  2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

ggplot(diamonds, aes(color, fill = cut)) +
  geom_bar(position = "dodge")